import java.util.Random;
import java.util.Arrays;
import java.lang.Math;

public class KStest
{
    public static void main(String[] args)
    {
	KStest ks = new KStest();
	ks.run();
    }

    public void run()
    {
	int n = 10000;
	int datarepeats = 1;
	int repeats = 1;
	Random r = new Random();
	double[] frac = {0.05, 0.1, 0.15, 0.20, 0.25};
	int fracs = frac.length;
	double[][] relErrPos = new double[3][fracs];
	double[][] relErrNeg = new double[3][fracs];

	for (int i = 0; i < 5; ++i)
	{
	    relErrPos[0][i] = 0.0; relErrPos[1][i] = 0.0; relErrPos[2][i] = 0.0; 
	    relErrNeg[0][i] = 0.0; relErrNeg[1][i] = 0.0; relErrNeg[2][i] = 0.0;
	}

	for (int k = 0; k < datarepeats; ++k)
	{
	    // initialize random data drawn from normal distribution N(0, 1)
	    double[] data = new double[n];
	    for (int i = 0; i < n; i++)
		data[i] = 2 * r.nextGaussian();   // changed stdev to 2
	    
	    // created sorted version of the data
	    double[] dataSorted = new double[n];
	    for (int i = 0; i < n; i++)
		dataSorted[i] = data[i];
	    Arrays.sort(dataSorted);
	    
	    // compute the real KS statistic
	    double KposReal = 0.0, KnegReal = 0.0;
	    for (int i = 0; i < n; i++)
	    {
		KposReal = Math.max(KposReal, (i + 1.0)/n - cdfNormal(dataSorted[i], 0, 1));
		KnegReal = Math.max(KnegReal, cdfNormal(dataSorted[i], 0, 1) - ((double)i)/n);
	    }
	    System.out.println("Real:\nKpos = " + KposReal + "\nKneg = " + KnegReal);
	    
	    int count = 0;
	    for (int l = 0; l < fracs; ++l)
	    {
		double f = frac[l];
		System.out.println("\nf = " + f);
		
		// test the GK-data structure
		double avgKpos = 0.0, avgKneg = 0.0;
		for (int j = 0; j < repeats; ++j)
		{
		    // insert data into GK data structure
		    GreenwaldKhannaQuantileSummary s = new GreenwaldKhannaQuantileSummary((int)(f * n));
		    for (int i = 0; i < n; ++i)
			s.insert(data[i]);
		    
		    // compute the estimated KS statistic from GK sketch
		    double Kpos = 0.0, Kneg = 0.0;
		    for (int i = 0; i < n; i++)
		    {
			Kpos = Math.max(Kpos, (i + 1.0)/n - cdfNormal(s.getQuantile(((double)i)/n), 0, 1));
			Kneg = Math.max(Kneg, cdfNormal(s.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);
		    }

		    avgKpos += Kpos/repeats;
		    avgKneg += Kneg/repeats;
		}
		System.out.print("GK: Kpos = " + avgKpos + ", Kneg = " + avgKneg);
		System.out.println("; " + relError(avgKpos, KposReal) + ", " + relError(avgKneg, KnegReal));
		relErrPos[0][count] += relError(avgKpos, KposReal)/datarepeats;
		relErrNeg[0][count] += relError(avgKneg, KnegReal)/datarepeats;


		// test the qDiqest data structure
		avgKpos = 0.0; avgKneg = 0.0;
		for (int j = 0; j < repeats; ++j)
		{
			// insert data into reservoir sampling data structure
			QDigest s = new QDigest (f * n/3.0);  // compression factor = memory/3
			for (int i = 0; i < n; ++i)
				s.insert(data[i]);
			
			// compute the estimated KS statistic from reservoir sketch
			double Kpos = 0.0, Kneg = 0.0;
			for (int i = 0; i < n; i++)
			{
				Kpos = Math.max(Kpos, (i + 1.0)/n - cdfNormal(s.getQuantile(((double)i)/n), 0, 1));
				Kneg = Math.max(Kneg, cdfNormal(s.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);
			}
				
			avgKpos += Kpos/repeats;
			avgKneg += Kneg/repeats;
		}
		System.out.print("QDigest: Kpos = " + avgKpos + ", Kneg = " + avgKneg);
		System.out.println("; " + relError(avgKpos, KposReal) + ", " + relError(avgKneg, KnegReal));
		relErrPos[2][count] += relError(avgKpos, KposReal)/datarepeats;                 
		relErrNeg[2][count] += relError(avgKneg, KnegReal)/datarepeats;
			
			
		// test the reservoir sampling data structure
		avgKpos = 0.0; avgKneg = 0.0;
		for (int j = 0; j < repeats; ++j)
		{
		    // insert data into reservoir sampling data structure
		    ReservoirSampler s = new ReservoirSampler ((int)(f * n));
		    for (int i = 0; i < n; ++i)
			s.insert(data[i]);
		    
		    // compute the estimated KS statistic from reservoir sketch
		    double Kpos = 0.0, Kneg = 0.0;
		    for (int i = 0; i < n; i++)
		    {
			    Kpos = Math.max(Kpos, (i + 1.0)/n - cdfNormal(s.getQuantile(((double)i)/n), 0, 1));
			    Kneg = Math.max(Kneg, cdfNormal(s.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);
		    }
		    
		    avgKpos += Kpos/repeats;
		    avgKneg += Kneg/repeats;
		}
		System.out.print("Sample: Kpos = " + avgKpos + ", Kneg = " + avgKneg);
		System.out.println("; " + relError(avgKpos, KposReal) + ", " + relError(avgKneg, KnegReal));
		relErrPos[1][count] += relError(avgKpos, KposReal)/datarepeats;                 
		relErrNeg[1][count] += relError(avgKneg, KnegReal)/datarepeats;
	    
		count++;
	    }
	    
	}

	for (int i = 0; i < fracs; ++i)
	{
	    System.out.println("f = " + frac[i]);
	    System.out.println("GK: " + relErrPos[0][i] + ", " + relErrNeg[0][i]);
	    System.out.println("QDigest: " + relErrPos[2][i] + ", " + relErrNeg[2][i]);
	    System.out.println("Sample: " + relErrPos[1][i] + ", " + relErrNeg[1][i]);
	}

    }

    // computes relative error of est from real
    public static double relError(double est, double real)
    {
	return Math.abs(est - real)/real;
    }
    
    // the cdf of N(mu, sigma) on value x
    public static double cdfNormal(double x, double mu, double sigma)
    {
	return 0.5 + 0.5 * erf((x - mu)/(sigma * Math.sqrt(2.0)));
    }

    // erf function taken from:
    // http://introcs.cs.princeton.edu/java/21function/ErrorFunction.java.html
    // fractional error in math formula less than 1.2 * 10 ^ -7.
    // although subject to catastrophic cancellation when z in very close to 0
    // from Chebyshev fitting formula for erf(z) from Numerical Recipes, 6.2
    public static double erf(double z) {
        double t = 1.0 / (1.0 + 0.5 * Math.abs(z));
        // use Horner's method
        double ans = 1 - t * Math.exp( -z*z   -   1.26551223 +
				       t * ( 1.00002368 +
					     t * ( 0.37409196 + 
						   t * ( 0.09678418 + 
							 t * (-0.18628806 + 
							      t * ( 0.27886807 + 
								    t * (-1.13520398 + 
									 t * ( 1.48851587 + 
									       t * (-0.82215223 + 
										    t * ( 0.17087277))))))))));
        if (z >= 0) return  ans;
        else        return -ans;
    }


}